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Abstract 

In this paper we consider the variabihty of the luminosity of a compact ob- 
ject (CO) powered by the accretion of an extremely inhomogeneous ("clumpy") 
stream of matter. The accretion of a single clump results in an X-ray flare: we ad- 
opt a simple model for the response of the CO to the accretion of a single clump, 
and derive a stochastic differential equation (SDE) for the accretion powered lu- 
minosity L{t). We put the SDE in the equivalent form of an equation for the 
flares' luminosity distribution (FLD), and discuss its solution in the stationary 
case. 

As study, we apply our formalism to the analysis of the FLDs of 

Super-Giant Fast X-ray Transients (SFXTs), a peculiar sub-class of High Mass 
X-ray Binary Systems (HMXBs). We compare our theoretical FLDs to the 
distributions observed in the SFXTs IGR J16479 - 4514, IGR J17544 - 2619 
and XTE J1739 - 302. 

Despite its simplicity, our model fairly agrees with the observed distributions, 
and allows to predict some properties of the stellar wind. Finally, we discuss how 
our model may explain the difference between the broad FLD of SFXTs and the 
much narrower distribution of persistent HMXBs. 

Subject headings: methods: analytical — methods: numerical — methods: 
miscellaneous — X-rays: binaries — X-rays: individual (IGR J16479 — 4514, 
IGR J17544 - 2619, XTE J1739 - 302 ) 
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1. Introduction 

In several astrophysical contexts the mass accretion onto a compact object (CO, 
a black hole, a neutron star a or white dwarf) cannot be considered as a continuous 
process. The accretion may be discrete: the mass does not flow as a continuous stream, 
but consists of a volley of clumps whose accretion on the compact object results in a 
flaring activity observed in the X-ray/7 I'ciy spectral window. Inhomogeneous flows occur 
in the Polar (a.k.a. AM Herculis) cataclysmic variables (e.g. Warner 1995), or in the 
Supergiant Fast X-ray Transients (Sguera et al. 2006, Negueruela et al. 2006; in't Zand 
2005). Understanding the properties of clumpy accretion on the super-massive black holes 
hosted at the core of galaxy clusters is also relevant in the context of the "cold feedback 
model" (e.g. Pizzolato and Soker 2005). 

The luminosity L{t) powered by the accretion on a compact object results from the 
interplay of two factors: the properties of the accreting stream and the response of the 
CO to their arrival. The response process may be quite complex, but it is essentially 
deterministic. On the other hand, if the accretion occurs from a population of clumps with 
random masses arriving at random times, the accretion driving term is a stochastic process. 
The accretion-powered luminosity L{t) is an irregular function of time, and it requires an 
adequate mathematical tool to be dealt with. This is provided by the theory of stochastic 
differential equations (SDEs), that replaces the rules of ordinary calculus to treat highly 
irregular functions (see e.g. Oksendal 2003; Gardiner 2009). 

In this paper we consider the variability of the X/7 ray luminosity of a CO powered 
by the accretion of a family of clumps with randomly distributed masses. We shall derive 
(Sec. 2) a simple SDE for the luminosity L(t), and will illustrate some of the subtleties 
involved in handling SDEs. A possible way to use this SDE is to compute a sample path 
of L{t), to be compared to an observed light curve. An equivalent approach (presented in 
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Sec. 3) consists in associating to the SDE for L{t) an integro-differential equation (known as 
"generalised Fokker-Planck equation" , GFPE) for the probabihty function p{L, t) that the 
source has luminosity L at the time t. We shall discuss the properties of the distribution 
p, and how it mirrors the features of the accreting stream. As a case study we apply our 
model to the super-giant Fast X-ray Transients (SFXTs), an interesting class of High 
Mass X-ray binaries composed by a neutron star (or a black hole) accreting an extremely 
inhomogeneous wind blown by its massive companion. Our analysis of SFXTs is presented 
in Sec. 4. We shall discuss our findings and in particular how our model helps to explain 
the difference of the FLDs of SFXTs from those of the much more common persistent X-ray 
sources. We shall summarise in Sec. 5. 

2. A simple stochastic accretion model 

The analysis of the long term variability of the X-ray light curve generated by the 
inhomogeneous accretion requires some modelling of the response of the CO the the 
accretion process of a single clump. To this purpose it is helpful to imagine the CO as a 
"black box" that emits X-rays in response to the accretion of mass from the surrounding 
environment. If Mc is the mass capture rate and L is the luminosity produced by the CO, 
we postulate the linear response 

L{t) = I dt'W{v,t-t')M,{t'), (1) 

where the function W (dependent on the parameters p) describes the response of the CO. 

The analytical shape of W may be exceedingly complex, since it must embody a number 
of physical processes: tidal effects, interaction of the accretion flow with the magnetic field 
of the CO, geometry and plasma instabilities of the accretion stream, radiation effects, and 
so on. Since this paper aims to focus on the basic properties of the accretion process, we 
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steer clear of these details and demand for W just few very basic properties; i) W{t) is 
causal: it vanishes for t < 0, i.e. before the capture of a clump, and ii) it decays to zero 
for t ^ 0, i.e. long after the clump has been captured. The simplest function with these 
properties is 

{0 t<0 

This response gives the accretion powered luminosity of a unit mass on a CO of mass 
M, down to a "stopping radius" R^t- If the CO is able to accrete the flow down to its 
radius (or event horizon) /?, then i?st = R- In some cases, e.g. a fast spinning magnetised 
neutron star, the so called "propeller" effect (see e.g. Lipunov 1987; Illarionov and Sunyaev 
1975; Davies et al. 1979; Davies and Pringle 1981) prevents the stream to flow beyond the 
magnetosphere, and in this case R^t is the Alfven radius Ra- The time scale r represents 
the time scale taken by the CO to "process" the accretion stream. For example, if the 
stream has an appreciable angular momentum, it forms an accretion disc before falling on 
the CO, and the mass accretion rate on the CO decays as a power law on the viscous time 
scale 

7-7.5x10 sc. 10-11 Me yr-ij \lA M,J era) ' 

(Eq. 5.63 in Frank et al. 2002), where -Rout is the outer radius and a < 1 is the viscosity 
parameter (Shakura and Sunyaev 1973). 

It is instructive to consider Eqns. (1) and (2) in the simple case of the accretion a 
single clump of mass m, with Mc{t) = m6(t) (where 6 is the Dirac delta function). With 
our choice (2) of the response W, the luminosity is zero at t < 0, i.e. before the clumps 
arrives, and 

Rst T 



- 6 - 



for t > 0. The luminosity has a sharp peak, and then decays exponentially with the 
e-folding time r. Plugging the Ansatz (2) into Eq. (1) we find 

dL GM . ^ 

We introduce the mass accretion rate M accreted down to the stopping radius Rst 

L . M, (6) 



so Eq. (5) becomes 



r^ = M.-A^. (7) 

at 



Since the function Mc{t) is random, the function M{t) is highly irregular and certainly not 
differentiable, so Eq. (7) is to be interpreted as a stochastic differential equation (see e.g. 
0ksendal 2003), and not as an ordinary differential equation. 

Thus far we have provided a simple model for the "response" of the compact object to 
the accretion of a clump of matter. We now turn our attention to the driving stochastic 
term M^, which embodies both the arrival rate of the clumps and their mass distribution. 
Neglecting the finite size of the clumps, we model the clumps' arrival rate as a train of delta 
pulses 

n{t) 

Mait) = J2^kSit-tk), (8) 

k=l 

where is the mass of the clump accreted at the time tk, and n{t) is a Poisson counting 
process, described by the probability 

r[n{t) =n] = (9) 

where A is the clumps' arrival rate. The masses {rrik} are distributed according to the 
probability distribution function ip{m)] we assume that ip is stationary, and that the random 
variables {m^} and {t^} are uncorrelated. The statistical properties of are readily 



derived: 

(Mc) = A (m) (10a) 
(M,(t) Me(t')) = ^ (^^) ^(^ - ^0- (10b) 

The properties (10) show that the Poisson process (8) is a white noise with non-zero mean. 

From Eqns. (7), (8) and (10) it is possible to derive some properties of the random 
function M{t). They not only provide a useful check against our more elaborate approach 
developed in the next sections, but they also help to point out some subtle key properties 
of SDEs. 

Integrating Eq. (7), we find the expression 

TdM = 5m{t) - M{t)dt, (11) 

where 

/t+dt 
dt'Mcit'). (12) 

With the aid of Eqns. (10) we derive 

/t+dt 
dt'{M^{t')) = X (m) dt (13a) 

i>t+dt pt+dt 



{{dmY) = J dt' J dt"{M,{t')M,{t")) = X {m') dt. (13b) 

With the first of these expressions we find the mean of Eq. (11), i.e. the ordinary differential 
equation 

r^ = X{m)-{M) (14) 
(the order of the mean and the derivative may be exchanged, see e.g. Reif 1965). 
For t ^ r the mean mass accretion rate (M) relaxes on the stationary value 

(M) = A(m). (15) 
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The evaluation of the second moment is trickier, on account of the stochastic nature of the 
equation. We write 

T d{M^) = T\{M + d Mf -mA =t \2MdM + (dM)^] , (16) 

where in the expansion we must keep the second order term {dM)"^ for reason that will be 
clear soon. Plugging Eq. (11) into this expression we get 

Td{M^) = 2M {5m- M{t)dt^ +{5mf + 0{dt), (17) 

and taking the mean by using Eqns. (13) we find the ordinary differential equation 

r = 2\^mf - 2 (M^) + - (m^). (18) 

dt T 

For t ^ T the second moment (M^) relaxes on the stationary value 

{M') = X^mf + ^{m'), (19) 

Z T 



corresponding to the variance 



al = {M')-{Mr = ^{m'). (20) 



Some remarks are in order. First, the stochastic nature of the function Mc makes M{t) an 
extremely irregular function of time. On account of such irregular behaviour, none of the 
rules of ordinary calculus are apphcable to M{t). One consequence is that the second order 
differential (dM)"^ is proportional to dt, and therefore it cannot be neglected in the formal 
manipulations, as we have done in Eq. (16). In order to work with such irregular functions 
a new kind of differential calculus (known as "Ito calculus" ) has been developed as a corner 
stone of the theory of stochastic differential equations (SDEs, see e.g. Oksendal 2003). 
Some level of knowledge of SDEs is therefore essential to master the modelling of random 
processes like that studied in this paper. 
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The moments (M") exist only if the correspondent moments (m") of the clumps' 
distribution are defined. This is not always the case, e.g. when the distribution ip has a fat 
tail. In order to deal with such possibility, it is necessary to extend the simple approach 
presented in this section. We shall explain how in Sec. 2.1. 

The mean accretion rate (15) depends on the clumps' arrival rate A, but not from the 
relaxation time r, which is an intrinsic property of the accretor. The variance (20) features 
the relaxation time r with an inverse power: a long relaxation time r corresponds to a 
narrow dispersion of the observed values of M around the mean. Indeed, if the relaxation 
time r is long, several clumps may be accreted before the accretor is able to respond, and 
an observer cannot distinguish between several elementary accretion processes, which are 
perceived as a single accretion of a large clump. Clearly, this has the effect of reducing the 
dispersion of the observed distribution. 

2.1. Direct SDE vs Fokker-Planck approach 

In principle, Eq. (7) may be solved numerically. Several computing techniques have 
been developed to tackle the numerical solution of SDKs (see e.g. Kloeden and Platen 
1999). A sample of the population of the clumps' masses and arrival times are extracted 
from the distributions ip and V[n{t)\ (defined by Eq. (9)). A sample path of the process 
M{t) is computed numerically, and M{t) is converted to the accretion luminosity L{t) via 
Eq. (6). A sample path Li{t) may be very different from another path L2{t), since the 
clumps' masses and arrival times extracted from the distributions (f and V[n{t)] used to 
compute Li{t) numerically may be very different from those extracted to compute L2{t). 
For this reason, it is meaningless to compare a single sample path of L{t) with a real 
light curve. Any comparison of a stochastic model with the observations will involve some 
statistics on a sample of the computed paths and the light curves. 
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We shall adopt an equivalent approach to achieve the goal of comparing a theoretical 
stochastic light curve with a real light curve. It is possible to associate to Eq. (7) a new 
equation (called a "generalised Fokker-Planck equation" , GFPE) for the probability p{L, t) 
that the source has luminosity L at the time t. The distribution p is directly comparable to 
the histogram of the flares' luminosities observed in a real source over a long time span. 

The approach based on the luminosity distribution p{L,t) allows to deal quite 
straightforwardly with a complication we have overlooked thus far. The luminosity 
distribution is defined as a function of the X-ray luminosity L, not of the mass accretion 
rate M. The general relation between the luminosity distribution p{L, t) and the mass 
accretion rate distribution P{M, t) is 

dM 



p{L,t) 



P{M,t) (21) 



dL 

If all the mass can reach the surface of the neutron star (located at the radius R), then 
Rst = R- In this case the distributions p{L) and P{M) are simply proportional to each 
other: 

p{m) = -^P{M,t). (22) 

In some systems, however, the accreting stream cannot reach the surface: this is typically 
the case of a magnetised fast spinning neutron star in which the matter is captured by the 
star's gravitational field, but is prevented from reaching its surface by a "propeller" barrier. 
At a very basic level, the "propeller" effect may be described as follows. The matter inside 
the magnetosphere (located at the Alfven radius Ra) is forced to rotate with the same 
angular velocity as the neutron star. On the other hand, the angular velocity of the 
mass stream outside the magnetosphere (r > Ra) is approximately keplerian, u ~ ujk{t). If 

Wns > (^k{Ra) (23) 

no accretion is possible, and the matter is stopped at the Alfven radius. The radius Ra 
depends on the geometry of the accretion flow. If NS accretes a clumpy wind in a detached 
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binary system the accretion is approximately spherical, and there are two possible regimes. 
The NS crosses a wind blown by the companion star with a velocity f^. This wind is 
focused by the NS's own gravitational field within the capture radius 

= ^ = 3.74 X 10" c,n ( ^) ( . (24) 

vl \1AMqJ VlO^kms-i/ ^ ^ 

If the magnetic field of the NS is "weak", then the wind is first focused by the NS's 
gravitational field, and only then it is affected by the NS' magnetic force. In this case the 
correct expression of the Alfven radius is (e.g. Lipunov 1987) 



= { ^^^^ ] ii Ra< Rg, (25) 



where /i is the magnetic moment of the NS. On the other hand, if the magnetic field is 
strong, then the wind feels NS' magnetic field before its gravity, and in this regime 



w 



Ra={ "^^/^ ) if Rc<Ra. (26) 



The critical mass accretion rate at which the propeller barrier sets in can be derived 
from Eq. (23) with Ra given by either Eq. (25) or Eq. (26). If M < M^^ the ram pressure 
exerted by the fiow is unable to overcome the magnetic pressure at the magnetospheric 
radius, and the fiow is stopped at the Alfven radius Ra- The stopping radius is then 

R M > 

Rst= { ■ (27) 

Ra M<M^ 

Plugging this stepwise function Rst{M) into Eq. (21) we find 

|^P(M,t) L<L, 

Li<L<L^, (28) 



p{L,t) 



GM 
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where Li = GMM^/Ra and = GMM^/R. 

For a typical parameters of a NS (/i ^ 10^*^ G) and of a stellar wind (f^ ~ 10^ kms~^) 
the luminosity Li is below ~ 10^^ ergs~^, which is barely observable. Therefore, we neglect 
the tail below Li and approximate Eq. (28) with 

{0 M <M^ 

. (29) 
^P(M,t) M>M^ 

In presence of an accretion barrier the luminosity distribution p{L) has a sharp lower cut-off 
at the luminosity L^^. 

Whether an accretion barrier is present or not, p must be computed from the mass 
accretion distribution P{M,t). If no accretion barrier is present, p is given by P{M,t) and 
Eq (22), else it is computed from P{M,t) and Eq. (29). 

In the next section we derive and discuss the GFPE for P{M,t) associated to the 
stochastic differential equation (7). 



3. The generalized Fokker— Planck equation 

The method to derive the GFPE from a stochastic differential equation is described by 
Denisov et al. (2009). We leave all the technical details of the derivation in Appendix A, 
and here we reproduce only the results derived there. The GFPE associated to Eq. (7) 
reads 

dP d{MP) 



- pP{M,t) + p dm^{m) p{m -m/T,t^ , (30a) 



dt dM 
where 

P = \t (30b) 
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is the accretion parameter. The solution of this equation is unique once we impose the 
normahsation 

POO 

/ dMP{M,t) = l (31) 
Jo 

and a suitable the initial condition. 



A dimensional analysis (e.g. Barenblatt 1996) helps to write Eq. (30a) in non- 
dimensional form. Suppose that ip depends on the characteristic mass mg as the only 
dimensional parameter. The only dimensional parameters in the distribution P are then 
M, t, T and mo. Note that the clumps' arrival rate A appears in the equation only in the 
dimensionless combination p = A r, and therefore it does not appear as a separate variable. 
The only way to combine the governing dimensional variables into dimensionless variables 
(from which P must depend) are 

x = Mr/mo s = t/T. (32) 

so 

p^pfMr \ 
J 

(we do not indicate the dependence on possible other dimensionless parameters). Plugging 
this expression into Eq. (30a) we retrieve 

dP dP ^ 

x — + {l-p)P + p / dy^{x-y) P{y,s), (34) 



ds dx jQ 

where P = Pttiq/t and (p = (pruQ are the scaled distributions P and ip. In the following we 
shall refer to the dimensionless Eq. (34), omitting the tildes on the scaled quantities. 



The problem (30) depends on some initial condition. After a long time t ^ t the 
compact star has gathered a large sample of clumps, and the probability P is expected to 
relax on an equilibrium, time independent configuration. In the remainder of this paper we 
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will shall assume that that such an equilibrium has been achieved, and we will work with 
the stationary distribution, solution of 

dP r 

X — = {p-l)P-p / dy ip{x - y) P (y) , (35a) 
ax Jq 

This is a Volterra integro-differential equation of the second kind, whose solution is unique 
once we impose the normalisation 

POD 

/ dxP{x) = l. (35b) 
Jo 

The solution to the problem (35) is particularly convenient via a Fourier (or Laplace) 
transform. This not only allows a numerical solution, but also yields some exact results on 
the moments of the distribution. 



3.1. Solutions to the generalised Fokker— Planck Equation 

In this section we derive a solution to the generalised Fokker-Planck equation. 
Although an exact solution does not exists for a general form of if, the Fourier transform of 
P is analytical. This allows to compute all the moments of P as function of the moments 
of (f (when they exist), and it is makes possible to compute numerically P{x) by evaluating 
a Fourier integral. 

The logarithm of the Fourier transform of Eq. (35a) reads 

lnP, = -p ^(l-V'.'), (36) 

where ipk is the Fourier transform of v?, and where we have introduced the normalisation (35b) 
as lnPfc=o = 0. Plugging Lpk into the last equation. 

In Pfc = p / dy(p{y)Em{iky), (37) 

^0 
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being 

Ein(C) = / ds (38) 

Jo ^ 

the entire exponential integral (see e.g. Abramowitz and Stegun 1972). The inverse Fourier 
transform 

f dk 

P{x) = / — ex-p[-ikx + In , (39) 

is rather involved, and cannot be further simplified for a general ip. Yet, it is possible to 
derive some exact results. The function InP^ is the cumulant distribution function (CDF) 
of P{x). We Taylor-expand the entire exponential integral in Eq. (37) and compare the 
resulting expression with the definition of CDF 

r=l 



lnP,^J2-^^^^ (40) 



to obtain the cumulants 



or, in dimensional units. 



= p {x')It, (41) 



((M.)> ^ ^ (42) 

From the cumulants it is easy to retrieve the first moments of the distribution P: in 
dimensional variables 



al^ = — {m') (43b) 



(M) = A (m) (43a 
A 

where (M), o"^^, 7^ and 72 are respectively mean, variance, skewness and excess kurtosis of 
P (see e.g. Sec. 14.1 of Press et al. 2007 for the definitions). 
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Mean and variance coincide with the values worked out in Sec. 2 from our simple 
analysis of Eq. (7). The higher order moments provide further qualitative details on the 
shape of the accretion distribution function P. The skewness 71 is always positive, meaning 
that P has a long tail for large x. The kurtosis 72 is always positive as well, i.e. P{x) is 
leptokurtic: it has a sharp peak and a heavy tail for large x. If the terms (m") are finite, 
the standardized moment of order n oi P scales as p^~"/^. 

The computation of the Fourier integral (39) cannot be done analytically for any ^p, but 
it can be tackled numerically. This will be done later for some choices of the distribution ip. 

If ip{m) is bounded for m — )■ 0, the asymptotic behaviour of P for x ^ 1 is independent 
of (p. For small x the integral in Eq. (35a) is negligible, and 

P{x) ~ xP-^ for a; < 1. (44) 

For small values of x, P is suppressed if p ^ 1, and diverges for p < 1. This behaviour has 
a clear physical explanation. For simplicity, we assume that there is no accretion barrier, so 
P may be read as a luminosity distribution. Consider the accretion of two small clumps, the 
second one hitting the system a time after the first one. The inequality p < 1 means 
that the relaxation time r is shorter than A~^: the compact object has time to process the 
first clump before the arrival of the second one, and their accretion triggers two successive 
flares. The accretion of each clump is "recorded" in the low luminosity tail, and therefore 
the value of P may be high. On the other hand, if p ^ 1, the compact object cannot 
process the first clump before the arrival of the second one. The accretion of two closely 
separated clumps cannot trigger two different low luminosity flares. Such accretion episodes 
are "recorded" as the accretion of a singe larger clump, and the low luminosity tail of P is 
depressed by this effect that for convenience will be referred to as "p suppression" . 
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3.2. A Dirac delta mass clumps' distribution 

Before applying the present model to a concrete astrophysical problem, it is useful to 
sketch the properties of P in the relatively simple case in which all the clumps have the 
same mass. The appropriate mass distribution is then (in dimensionless units) 

(/.(x) = 5(x-l), (45) 

where 6 is the Dirac delta function. The generalised Fokker-Planck equation (35a) reduces 
to the delay differential equation 

x'^ = {p-l)P -'dix-l) Pix-l), (46) 
dx 

where ^{u) is the Heaviside step function, -i? = for n < and ^9 = 1 for n > 0. The 
solution of Eq. (46) can be easily tackled numerically. Some plots of P{x) are shown in 
Fig. 1 for several values of the accretion parameter p. 

03 
OA 
0.3 
0.2 
0.1 
O.Q 

0.01 0.10 1.00 10.00 

Fig. 1. — The accretion PDF for the delta clumps' distribution (45) for several values of the 
accretion parameter p. 
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The mean of P is (x) = p and the standard deviation is o" = a/ p/2. On account of 
the lack of massive clumps in the distribution (45) the distribution P does not display a 
significant tail. If p > 1, P is also suppressed for x < 1. These effects make the distribution 
P quite narrow around its mean. The limit distribution solution of (46) will be useful later 
as a benchmark for more complex cases. 



V^ix) = = exp 



(47) 



3.3. A log-normal clumps' mass distribution 

We now explore in some detail the solution of the GFPE (30a) with a more general 
distribution of the clumps' masses. To our purposes, it is particularly convenient to adopt 
the log-normal distribution 

In^ X 

defined for any a; > 0. This choice has several advantages. First, the log-normal is extremely 
flexible: by tuning the shape parameter a, it may mimic distributions as different as a Dirac 
delta (for cr ^ 1), or a power law (for cr > 1). The distribution (47) may be written (up to 
a constant factor) if ~ x~^, where ( = 1 + lnx/2cr^. The exponent C is a slowly varying 
function of the range x, and for x <^ e'^^/^ the log- normal (f closely resembles a power law 
with index C ^ 1- All the moments of a log-normal are defined, being 

(x'^) =exp(nVV2). (48) 

Eqns. (43) and (48) give the moments of the distribution P 

(M) = Amoe"'/^ (49a) 
al, = ^ml e-^ (49b) 

7i = g|e-^/^ (49c) 

2 



72 = e^" /p. (49d) 
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It is clear that the moments of P are defined for any order ra, but their magnitude rapidly 
increases with n ii a > 1, making P heavy tailed. 

The Fourier integral (39) cannot be computed exactly if </? is a log-normal, but it must 
be evaluated numerically. The numerical calculation of Fourier integrals poses a well known 
problem, widely discussed in the mathematical literature. A very convenient numerical 
method is a variant of the double exponential quadrature (see e.g. Mori and Sugihara 2001) 
suggested by Ooura and Mori (1999). 

In Fig. 2 we show some theoretical curves computed from a log-normal ip for several 
values of a and p. If p is fixed (see the top panel of Fig. 2), as a increases, also the 
distribution becomes increasingly broad, stretching its tail to high values of x. The peak of 
P moves to the right with increasing a, and its height lowers to preserve the normalisation. 
For moderately large values of cr, the interval over which P significantly differs from zero 
spans several orders of magnitude. 

For a fixed a, on the other hand, (bottom panel of Fig 2) for increasing values of p 
the distribution P is suppressed for low x. This general property of P has been already 
discussed at the end of Sec. 3.1. 

4. Application to Supergiant Fast X-ray Transients 

We apply the formalism developed in this paper to the analysis of the fiares' distribution 
observed in the Supergiant Fast X-ray Transients (SFXTs). We first outline the main 
properties of this class of objects, and then will analyse them in the framework of stochastic 
accretion. 



i 

0.0001 0.0010 0.0100 0.1000 1.0000 10.0000 100.0000 



Fig. 2. — Two samples of theoretical curves computed from the log-normal distribution (47). 
In the top panel p is fixed at p = 10 and the curves refer at different values of a. In the 
bottom panel a has been fixed at a = 1 and the curves refer to different values of p. 
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4.1. Presentation of SFXTs 

SFXTs are a class of rare transient X-ray sources (about ten members are known up 
to date), associated to OB supergiant stars (see Sidoli 2011 for a recent review), making 
them a sub-class of the HMXBs. 

The short and bright X-ray emission from several members of this class was discovered 
by the INTEGRAL satellite during monitoring observations of the Galactic plane at hard 
X-ray energies (Sguera et al. 2006, Negueruela et al. 2006). About a half of them are X-ray 
pulsars and the mass flow from the massive donor to the compact star (usually assumed to 
be a neutron star in all the members of the class) occurs via a stellar wind. 

SFXTs spend most of their time in a relatively low level brightness state (below 
Lx ~ 10^^ ergs"^; Sidoli et al. 2008), reaching Lx ~ 10^^ ergs~^ in quiescence, but show 
an occasional "outbursting" activity when the mean X-ray luminosity is much higher than 
usual for a few days, and is punctuated by a sequence of short (from a few minutes to a few 
hours long) flares peaking at X-ray luminosities of 10'^^ — 10'^'' ergs~^ (e.g. Romano et al. 
2007; Rampy et al. 2009). 

Their X-ray transient activity with such a wide dynamic range (from 2 to 5 orders of 
magnitude) is puzzling, since the components of these binary systems seem very similar to 
persistently accreting classical X-ray pulsars (e.g. Vela X-1). Two kind of explanations 
have been proposed, although none of them is able to completely explain the whole 
phenomenology of all the members of the class, to date. 

In the flrst type of models, the compact object accretes mass from an extremely 
inhomogeneous "clumpy" stellar wind (in't Zand 2005), which is believed to be present in 
massive OB stars (see e.g. the review of Puis et al. 2008). In a few cases a preferential plane 
for the outflowing clumpy wind is suggested (Sidoli et al. 2007, Drave et al. 2010) while 
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usually a spherically symmetric morphology is assumed (Negueruela et al. 2008). It is still 
unclear if the flow is able to form a (possibly transient) accretion disc (see e.g. Ducci et al. 
2010). 

In the second kind of models, the compact object is a slowly spinning magnetic neutron 
star undergoing transition between the accretion regime and the onset of a centrifugal or a 
magnetic barrier (Grebenev and Sunyaev 2007, Bozzo et al. 2008). This scenario involves 
the interplay of a host of physical processes, including several kinds of hydrodynamical 
instabilities the details of which are beyond the simplified sketch outlined in this paper. 

A sample of SFXTs has been extensively observed by the X-ray satellite Swift 
(Romano et al. 2010). Fig. 3 shows the histograms of the flares' luminosity distributions 
for three SFXTs that we computed from the data gathered in a two years' campaign with 
Swift (Romano et al. 2010). 

These histograms are close relatives to the flares' luminosity distribution FLD 
introduced in Sec. 1. We derived the X-ray luminosity histograms from the distributions 
of count rates originally reported in literature by Romano et al. (2010). Conversion factors 
from Swift XRT count rates to X-ray luminosities (1 — 10 keV) have been calculated 
assuming a power law spectrum, with the distances and the spectral parameters appropriate 
for each individual SFXT, as reported in Sidoli et al. (2008). 

4.2. The model 

In this section we interpret the observed luminosity distributions plotted in Fig. 3 in 
the light of the stochastic set up in the previous sections. Some of the questions we address 
are the following: 

• The SFXTs' luminosity distribution drops quite abruptly below Lx — 10^^ — 
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Fig. 3. — Distribution of the X-ray luminosity of three SFXTs we calculated from the X-ray 
count rate distribution observed with Swift, during a two years' campaign. Data have been 
taken from Romano et al. (2010). 

Iqsa gj-gg-i Provided that this is not a selection effect in our observations, what is 
the origin of this feature? 

• What properties of the accretion stream and accretion process can we infer by 
comparing our model with the observations? How do the properties of the stellar 
wind we deduce compare with the existing literature? 

• Can we explain the difference between the wide dynamic range of luminosities 
observed in SFXTs and the much narrower one seen in persistent sources, like e.g. 
Vela X-1? 

Our comparison of the model with the observed flares' luminosity distributions (FLDs) 
is not intended to be a "best fitting" procedure, for several reasons. 

First of all, our model (2) for the response of the accretor may be a bit oversimplified. 
The theoretical luminosity distribution P computed from it, therefore, is not expected to 
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provide an accurate model of the SFXT flares' luminosity distribution observed in the real 
world ^. Second, only the FLDs of three sources are available for our analysis, probably too 
few to draw definite conclusions. Third, not all the flares may have been spotted during 
the observation campaign, and the FLDs we work with might not describe a fair sample 
of all the flares. Finally, the calculation of a model for given values of p and a is quite 
time consuming, and a quantitative determination of those parameters via a best fitting 
procedure is not viable at the moment. 

When we compare an observed FLD with its theoretical model, we aim to get a hint at 
the properties of the accretion process (described by the parameter p) and the properties 
of the family of clumps (described by the log-normal distribution (47)). A more physically 
motivated response is clearly necessary in order to make more quantitatively definite 
statements. 

Despite its all too obvious limitations, we believe that our model has its own value. 
Its qualitative predictions are expected to hold even for a more sophisticated model; in 
addition it points out the fundamental properties of stochastic accretion that can be easily 
shaded by rather opaque mathematical complications that we leave to a forthcoming paper. 
That said, we are now ready to compare out theoretical predictions with the observed 
flares' distributions in the three SFXTs IGR J16479 - 4514, IGR J17544 - 2619 and 
XTE J1739 - 302. 

^Maybe the response (2) is not always thus bad. The SFXT IGR J16479 - 4514, for 
instance, is known to exhibit flares that after a fast rise decay exponentially (Sguera et al. 
2005; Ducci et al. 2010). Sguera et al. (2005) could fit the light curve of a flare with sufficient 
statistics with an exponential law with an e-folding time r = 15 ± 9 min. 
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We aim at determining the best fitting parameters of our model from the flares' 
distributions of a given SFXT. In general, the free parameters of the model are the accretion 
parameter p, the log-normal shape factor a, the ratio m^jr and the cut-off mass accretion 
rate due to the action of an accretion barrier (see Eq (29)). In order to keep things as 
simple as possible, in the following analysis we assume that no accretion barrier is present 
(i.e. = 0), and the the model is completely defined by the three parameters p, a and 
mo/r. 

What is the best way to determine these parameters from the data of an observation 
campaign of a given SFXT? As we shall see, our choice is rather limited. 

One might think, for instance, to determine the parameters by comparing the mean, 
variance, skewness and excess kurtosis computed from our model (see Eq. (43)) with the 
observed mean, variance, skewness and excess kurtosis, provided that we are able to convert 
the mass accretion rate M to a X-ray luminosity. If there is no accretion barrier, the 
modelled luminosity distribution is simply "piV) = P(x)/L*, where P{x) is the dimensionless 
distribution computed from the generalized Fokker-Planck equation, and 

G niQ 

L* = p ^ , (50) 

where M^s ~ 1-4 Mq and R^s ~ 10 km are the typical mass and radius of a neutron star. 
A direct comparison of the moments of the theoretical and the computed distributions, 
however, is not viable. The use of high order momenta such as the skewness or the kurtosis 
of a heavy tailed distribution (like our P{x)) is not recommended for any statistical analysis 
(e.g. Press et al. 2007). The problem is that these moments are not robust, i.e. they are too 
sensitive to the outliers. Any small change in the number of the observed high luminosity 
flares would alter the values of these moments considerably, making them of little use to 
any statistical purpose. 
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Since we cannot use the moments, we are forced to determine the parameters by 
comparing the full shapes of the observed and the theoretical distributions. The data from 
Romano et al. (2010) available to our analysis are already binned in luminosity intervals: 
we known the fraction pj^°^^'^ of the flares observed in the luminosity interval but 
we don't have access to the luminosity of each single flare. This prevents us from using any 
variant of the Kolmogorov-Smirnov test to compare the theoretical and the observed flares' 
distributions. 

For these reasons, we must resort to a least squares method to determine the 
parameters p, a and rrio/T. Suppose for the moment that p and a are known, and the only 
parameter to be determined is ttiq/t. For a (tentative) value of mo/r we compute the scale 
luminosity (50) and the expected fraction pj^^^^^ of flares within the luminosity interval 
[Lj,Lj+i]. We compute the chi-squared 



where we assume the errors on the observed fractions ~ Pj . Finally, we find the value 
of mo/r that minimises X^. A minimisation procedure to determine all the parameters p, a 
and niQ/r would be very time consuming, and therefore is not viable at the moment. For 
this reason we have built a grid of values of p and a, and for each couple (p, a) we have 
computed the correspondent P{x), and the value of mo/r that minimises the for that 
P. Finally, we chose as our best estimates for p, a and mo/r the values that return the 
smallest X^. On account of the (expected) poor quality of the fits, we do not report the 
values of the X"^ statistics, but we only discuss here below our results. 

Fig. 4 shows the flares' distributions for our three SFXTs superposed to their "best 
fitting" models. On account of the asymptotic behaviour (44) of the theoretical luminosity 
distribution, a glance to Fig. 3 hints that in all cases it must be p > 1 to match the observed 




(51) 



.2 
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FLDs. In addition, the observed high dynamical ranges of the flares' luminosities require a 
heavy tailed clumps' mass distribution, i.e. a relatively high log-normal shape parameter a. 

IGR J17544 — 2619 is the source with the lowest dynamic range, since it displayed 
flares ranging in the luminosity interval between 3.9 x 10^^ ergs~^ and 2.7 x 10'^^ ergs~^. 
There is no apparent low luminosity cut-off (which might be interpreted as the effect of an 
accretion barrier). The FLD of this source may be reproduced quite nicely by a model with 
p = 8, a = 4 and tjiq/t = 1.0 x 10^^ g s"^. 

The source IGR J16479 — 4514 has the highest dynamic range of all the considered 
sources, with flares' luminosities ranging from 1.4 x 10^^ ergs"^ to 3.5 x 10^'' ergs"^. In this 
case the best value of the log-normal shape is a = 5 (slightly larger than the value found for 
IGR J17544 — 2619). The best estimate of the accretion parameter is p = 10, and even in 
this case it is not probably necessary to invoke an accretion barrier to explain the left tail 
of the observed FLD. The last parameter is mo/r = 5.5 x 10^^ g s~^. 

The flares of XTE J1739 — 302 show luminosities between 3.3 x 10^'^ ergs~^ and 
5.3 X 10^^ ergs"^ Its steep "wall" that cuts off S{L) below L = 2.1 x lO^^ ergs"^ is difficult 
to reproduce. Following our discussion of the "p suppression" at the end of Sec. 3.1, one 
might expect that a very high value of p is necessary to cut off abruptly the low tail of the 
luminosity distribution. Indeed, our best-fitting model for this source gives p = 40, a = 8 
and mo/r = 1.0 x 10'' g s~^. We note however that the space of parameters characterised 
by high values of p and a is quite difficult to explore with our means. For this reason, and 
also on account of the relatively poor quality of the "best fit" shown at the bottom panel 
of Fig. 4, we cannot be sure that the p suppression is winning over the "propeller barrier" 
process discussed in Sec. 2.1. 




Fig. 4.— The observed FLDs of the SFXTs IGR J17544-2619 (top panel), IGR J16479- 
4514 (middle panel) and XTE J1739 — 302 (bottom panel) superposed to their "best fitting" 
theoretical curves. The closest agreement between the observed and the computed FLD are 
found for p = 8, (T = 4 and mo/r = 1.0 x 10^^ g g-i_ j^fo^ IGR J17544 - 2619), p = 10, 
(7 = 5 and mo/r = 5.5 x lO^^ g s'^ (for IGR J16479 - 4514) and p = 40, a = 8 and 
mo/r = 1.0 X 10^ g s^^ (for XTE J1739 - 302). 
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To resume, we may reproduce at least the gross properties of the sources' FLDs, with 
moderate values of the accretion parameter p and a relatively high shape parameters a of 
the clumps' log-normal. At least two of our sources do not require any accretion barrier to 
explain their low luminosity tail, that may be justified by the "p suppression". The source 
XTE J1739 — 302 is more problematic, and here a propeller barrier is not excluded. 

Our model has three degrees of freedom: mo/r, a and the accretion parameter p = \ t. 
The mass arrival rate A is degenerate since it only appears in a product with r in the 
accretion parameter p, and so it cannot be determined directly from a comparison of our 
model with the observations. We may however get a hint from the analysis of the flares' 
light curves and from the geometry of the system. 

In their analysis of the light curve of IGR J16479 — 4514 Sguera et al. (2005) found 
that the flares lasted from ~ 30 min to ~ 3 h. In one case, they were able to analyse 
a single light curve that after a fast rise, decayed exponentially with an e-folding time 
r = 15 ± 9 min. Exponentially decaying light curves are not uncommon for the flares 
occurring in this source (Ducci et al. 2010). If the decay time r ~ 10"^ s may be taken as 
typical, then for p ~ 10 we estimate the time interval between the successive accretion of 
two clumps ~ 10^ s. If the NS orbits a massive star of M ~ 20 Mq with a period 
Porb — 3d (typical of IGR J16479 — 4514), the distance covered by the NS between two 
successive encounters is of the order of few 10^ cm, a fraction of the stellar radius. From 
the values of mo/r and r ~ 10'^ s we find values of mo between ~ 10^° g and ~ 10^^ g for 
our three sources. On account of the large values of a, the mean masses are much larger 
than these (compare with Eq. (48)), and are of the order of 10^^ — 10^^ g. These values for 
the clumpy wind are in broad agreement with estimates obtained with other means (see 
e.g. Ducci et al. 2009). 
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We now present a simple argument showing how the parameters vtlqIt^ a and p may 
not be fully independent of each other. Consider a giant donor star with mass and 
radius i?* that blows a clumpy wind. Assuming that the mass carried by the wind is all in 
clumps, then the mean stellar mass outflow rate at the distance r from the star is 

= 4 7rr^nciWci (m), (52) 

where is the mean number of clumps per unit volume, and v^x is the mean clumps' 
velocity. As the neutron star orbits the giant star, it crosses this clumpy wind and encounter 
clumps with the mean rate 

\ = -K R% Worb n^h (53) 

where forb is the orbital velocity of the NS and Rg is the capture radius defined by Eq. (24). 
Eliminating from Eqns. (52) and (53) and postulating that the population of the masses 
is log-normally distributed with shape parameter a and median mass mo we find 

A= f^V M^'!^^-<rV2 (54) 
V 2 r y niQ Vc\ 

The rate A is a strongly decreasing function of the clumps' mass distribution shape a for 
the following simple reason. If a is large, the log- normal distribution is skewed to the right, 
and most of the mass blown by the wind is carried by large, heavy clumps. A given is 
carried by relatively few clumps, and their encounters with the orbiting neutron star are 
rare, i.e. A is small. Multiplying the last equation by r, we find 

p = f ^ ^ e-'^V2. (55) 

\2r J mo/T Vd 

This shows that in our model the parameters a, p and mo/r are not independent of each 
other: p is proportional to (mg/r)"^ and decreases exponentially with It is difficult to 

use this equation to make quantitative predictions about the magnitude of p as a function 
of the other parameters. Indeed, when made explicit, its right hand side is found to be very 
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sensitive to the values of rather uncertain quantities (e.g. the velocity of the stellar wind), 
and even a small change in their value may vary p by orders of magnitude. We may plug 
the values of A and p into Eqns. (49) to find 



where (L) and are respectively the mean accretion luminosity and the standard deviation 



an explanation for the difference between the ordinary, persistent X-ray sources occurring 
in the high mass binaries, and the SFXTs. Indeed, according to Eq. (56), even moderate 
values of a (say, a ^ 4 — 5) greatly amplify the dispersion around the mean (Jl/ (L). This 
means that the wide dynamical range of luminosities observed in SFXTs is but the effect 
of the accretion of a population of clumps with a wide mass spectrum (f{m). On the other 
hand, if a ^ 1, i.e. all the accreting clumps have very similar masses, there is no such 
strong amplification, and aL/{L) may be small. The FLD generated by the accretion of 
clumps with narrowly distributed masses should resemble the distribution P discussed in 
Sec. 3.2. 

According to this model, then, the difference between persistent HMXBs and SFXTs is 
entirely due to a substantial difference between the winds blown by the donor stars hosted 
in these systems. Since SFXTs are much rarer than the ordinary HMXBs, probably the 
occurrence of a clumpy stellar wind with a wide range of masses is also a rare event. 




(56) 



of 



L, if L oc M. The exponential factor at the right hand side of the Eq. (56) suggests 



5. Summary and conclusions 



In this paper we have set up a simple stochastic model for the accretion of a clumpy 
stream on a compact object. We have modelled the response of the compact object to the 
accretion with a simple exponential law (2) characterised by the relaxation time r. The 
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accretion stream is described as a train of pulses hitting the compact object with the mean 
rate A. The clumps have a mass distribution described by the function ip{m). With these 
ingredients we have written a stochastic differential equation (SDE) for the mass accretion 
rate M{t). We have also derived the generalised Fokker-Planck equation (GFPE) associated 
to the SDE. Its solution P{M, t) is the probability that the mass accretion rate is M at 
the time t. We have then restricted us to the analysis of the stationary solution P{M). 
In general P cannot be written in a closed form for any choice of ip, but it is nevertheless 
possible to express its moments as functions of the moments of ip (if they exist). In addition, 
the asymptotic behaviour of P{M) for small M is P{M) oc Af'^, for any (f{m) limited for 
small m. For p ^ 1, then, P{M) is negligible for M < ttiq/t. This "p suppression" of P 
occurs when the relaxation time r is long with respect to the clumps' arrival rate, and the 
compact object cannot respond promptly enough to the rapid succession of clumps' arrivals 
(Sec. 3.1). 

We have then studied in some detail the function P for two choices of v?: a delta 
function and a log-normal distribution. In the case of a log-normal (f, P depends on three 
parameters: the accretion parameter p ^ Xr, mo/r and a, where mo and a are the location 
and shape parameters of the log-normal. 

As a case study, we have applied our formalism (with a log-normal (/?) to the physics 
of Super Giant Fast X-ray Transients (SFXTs), a peculiar sub-class of High Mass X-ray 
binaries (HMXBs). 

We have compared our model to the flares' X-ray luminosity distributions of the 
SFXTs IGR J17544 - 2619, IGR J16479 - 4514 and XTE J1739 - 302 observed in a two 
years' campaign with Swift. 

We have found that, however simplified, our model may reproduce the features of the 
observed flares' distributions. The typical parameters giving the best agreement between 
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the model and the observations are p ~ 10 and a ~ 5, corresponding to a long tail in the 
mass distributions. The parameters r and A cannot be determined separately, but taking 
the value r ^ 10^ s suggested by some observations, we find ~ 10^ s, or an inter-clump 
average distance 10^ — 10^° cm. The average masses of the clumps are in the order of 
10^*^ — 10^^ g, in broad agreement with the estimates found in the existing literature (e.g. 
Ducci et al. 2009). In at least two of the SFXTs we analysed, we were able to reproduce 
the fiares' luminosity distributions without invoking an accretion barrier. 

The reason why SFXT show a much higher dynamic range than the ordinary HMXBs is 
still unclear: according to our model, the difference is entirely due to the different properties 
of the streams of matter accreting the compact objects hosted in those systems. HMXBs 
accrete from clumpy winds with a narrow mass distribution, while SFXTs from winds with 
an extremely wide clumps' mass distribution. In a recent paper Oskinova et al. (2012) 
pointed out that a stochastic component in the clumps' velocities is probably essential in 
shaping the light curves of SFXTs. In the present paper we have neglected the possible 
contributions to L{t) from a stochastic component of the velocity of the accreting stream, 
but it must clearly be included in any (more) realistic model. 



In this section we closely follow Denisov et al. (2009) to derive the generalized 
Fokker-Planck equation (30a) from the Langevin Equation (7). According to the Ito 
interpretation, the solution of Eq. (7) in the (short) time interval At is 



A. Derivation of the generalized Fokker— Planck Equation 



M{t + At) 



M{t)-At M{t)/T + 6M^{t)/T, 



(Al) 



where 




(A2) 
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being n{At) is the number of clumps accreted in the time interval At. 

In order to derive the generalised Fokker-Planck equation we introduce the probability 
p(AM, At) that the mass AM is accreted by the compact object in the interval At 

p{AM, At) = {5{AM - SM,{t))). (A3) 

From this definition and Eq. (A3), after some algebra it can be shown that 

p{AM, At) = Po(At) 6{AM) + W{AM, At) (A4a) 

where 



W{AM, At) = ^ Pk{At) / dmi ip{mi) / dm2 </?(m2) ■ ■ ■ / druk (p{mk) 8 ( AM - ^ 
k=i J J J \ 



k 

777,1 



(A4b) 

being P„(At) the Poisson probability that n clumps are accreted in the time interval At, 
and ip{m) is the mass distribution of the clumps. If At is small, we may expand p to first 
order in At and take 

p(AM, At) = (1 - AAt) 5{AM) + AAt Lp{AM) + C(At) (A5) 

The accretion probability P{M, t) is defined as 

P(M, t) = {6{M - M{t)). (A6) 

In order to proceed, we need a couple of formulae to express the averages of 
the functions {F{M(t))) and {F{M(t),6Mc{t))) in terms of the accretion probability 
distributions P{M,t) and p(AM, At). Since the random variable M(t) and 6Mc{t) are 
independent, then 

(F(M(t))) = [ dM P{M,t) F{M) (A7a) 



and 



{F{M{t),SMc{t))) = j dM j d{AM) P{M,t) p{AM,At) F(M,AM). (A7b) 

We now introduce the time Fourier transformation to derive the generahsed Fokker-Planck 
equation. From the definition it is clear that 



Pfc(t) = / dM P{M,t) e 



-ikM 



We calculate the increment SPk{t) = Pk{t + At) — Pk{t) from Eq. (Al) as 



(A8) 



(A9) 



Expressing M{t + At) with Eq. (Al) and retaining only the terms linear in At, 



SPkit) 



-ikM{t) 



-ikSMc{t)/T 



At 

T 



+ ik^{M{t) e-*'^^''(*A + 0{At) 



(AlO) 



Introducing the probability distributions P and p, 



5Pk{t) 



dM P{M,t)e-'^^ I d{AM)p{AM,At) 



+ik— I dMP{M,t)Me-'^^ +0{At) 



-ikSMc{t)/7 



+ (All) 
(A12) 



Plugging Eq. (A5) into this expression we obtain 



5Pk{t) 



dM P{M,t)e-'^^'XAt d{AM)^{AM) 



-ikSMc{t)/T 



At 



+ik— I dMP{M,t)M e-'^^^ + 0{At) 



Dividing by At and sending At — )■ 
dPk 



XPk Vk/r -XPk + i- I dMP{M, t)M e 



-ikM 



Finally, the inverse Fourier transform yields 

— = -XP+X dm<^(m) P (M -m/T,t) + ^{MP 

dt Jo J y ' ) T dM 



(A13) 
(A14) 



(A15) 



(A16) 
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which is the desired Fokker-Planck equation. 
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